Many-body physics in the radio frequency spectrum of lattice bosons 
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We calculate the radio-frequency spectrum of a trapped cloud of cold bosonic atoms in an optical 
lattice. Using random phase and local density approximations we produce both trap averaged and 
spatially resolved spectra, identifying simple features in the spectra that reveal information about 
both superfluidity and correlations. Our approach is exact in the deep Mott limit and in the deep 
superfluid when the hopping rates for the two internal spin states are equal. It contains final state 
interactions, obeys the Ward identities (and the associated conservation laws), and satisfies the f- 
sum rule. Motivated by earlier work by Sun, Lannert, and Vishveshwara [Phys. Rev. A 79, 043422 
(2009)], we also discuss the features which arise in a spin-dependent optical lattice. 
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I. INTRODUCTION 



Bosonic atoms in optical lattices, described by the 
Bose- Hubbard model p], display a non-trivial quan- 
tum phase transition between a superfluid and Mott in- 
sulator. The latter is an incompressible state with an 
integer number of atoms per site. In a trap the phase 
diagram is revealed by the spatial structure of the gas: 
one has concentric superfluid and insulating shells. This 
structure has been elegantly explored by radio frequency 
(RF) spectroscopy a technique which has also given 
insight into strongly interacting Fermi gases across the 
BEC-BCS crossover Here we use a Random Phase 
Approximation (RPA) that treats fluctuations around 
the strong coupling Gutzwiller mean field theory to ex- 
plore the radio-frequency spectrum of lattice bosons. 

We find two key results: (1) Our previous sum-rule 
based analysis Q of experiments at MIT Q stands up 
to more rigorous analysis: in the limit of small spectral 
shifts, the RPA calculation reduces to that simpler the- 
ory. (2) In a gas with more disparate initial and final 
state interactions (such as Cesium), the spectrum be- 
comes more complex, with a bimodal spectrum appear- 
ing even in a homogeneous gas. The bimodality reveals 
key features of the many-body state. For example, in the 
limit considered by Sun, Lannert, and Vishveshwara @, 
the spectral features are related to the nearest-neighbor 
phase coherence. In the Gutzwiller approximation, the 
phase coherence directly maps onto the condensate den- 
sity. In this paper we provide a physical picture of this 
result and explain how this bimodality can be observed 
in a spatially resolved experiment. 



A. RF Spectroscopy 

In RF spectroscopy, a radio wave is used to flip the 
hyperfine spin of an atom from \a) to \b). The rate of 
excitation reveals details about the many-body state be- 
cause the | a) and \b) atoms have slightly different in- 



teractions. Generically the interaction Hamiltonian is 
Hint = J2j U aa n a (n a - l)/2 + U bb n b (n b - l)/2 + U ab n a n b , 
with U aa T U ab 7^ Ubb, where n a is the number of a-state 
atoms on site j. In the simplest mean- field picture, the 
energy needed to flip an atom on site j from state a to 
state b is shifted by an energy 5uo = U bb n b + (U ab — U aa )n a . 
Applying this picture to an inhomogeneous gas suggests 
that the absorption spectrum reveals a histogram of the 
atomic density. Such a density probe is quite valuable: 
in addition to the aforementioned examples, it was the 
primary means of identifying Bose-Einstein condensation 
in atomic hydrogen Q. 
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FIG. 1: (Color online) Illustration of two types of RF-active 
excitations of the lattice superfluid near the Mott transition. 
Open (blue) circles are atoms in the \a) state, filled (red) 
circles are atoms in the \b) state, and the arrows indicate a 
delocalized particle while other particles are localized, (a) 
Illustrates the initial superfluid state, consisting of a dilute 
gas of atoms moving in a Mott background. Final states in 
(b) and (c), show the excitation of a core or delocalized atom. 

Recently Sun, Lannert, and Vishveshwara found a 
bimodal spectrum in a special limit of this problem, as 
did Ohashi, Kitaura, and Matsumoto [8] in a separate 
limit, calling into question this simple picture. We give 
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a simple physical interpretation of the bimodality. As il- 
lustrated in Fig. [TJ the superfluid state near the Mott in- 
sulator can be caricatured as a dilute gas of atoms/holes 
moving in a Mott background. An RF photon can ei- 
ther flip the spin of one of the core atoms, or flip the 
spin of one of the mobile atoms. The energy of these two 
excitations will be very different, implying that the RF 
spectrum should be bimodal. Through our RPA calcula- 
tion, we verify this feature, calculating the frequencies of 
the two peaks and their spectral weights. Interestingly, 
this calculation reveals that the two excitations in our 
cartoon model are strongly hybridized. 

We find that that for parameters relevant to experi- 
ments on 87 Rb, that the degree of bimodality is vanish- 
ingly small and our previous sum rule arguments [H ac- 
curately describe such experiments. On the other hand, 
there are opportunities to study other atoms (for exam- 
ple, Na, Cs, Yb) for which the bimodality may be more 
pronounced. Moreover, if the interactions or tunneling 
rates can be tuned via a spin-dependent lattice or a Fes- 
hbach resonance then this spectral feature will appear in 
a dramatic fashion. 

This bimodal spectrum, with one peak produced by the 
"Mott" component and another by the "superfluid" com- 
ponent, is reminiscent of the spectrum of a finite temper- 
ature Bose gas in the absence of a lattice. As described 
by Oktel and Levitov [9], in that situation one sees one 
peak from the condensate, and one from the incoherent 
thermal atoms. We would expect that at finite temper- 
ature our "Mott" peak continuously evolves into their 
"thermal" peak. 



II. BOSE-HUBBARD MODEL 
A. Model and RF spectra 

In the rf spectra experiments we consider, initially 
all atoms are in the a-internal state and the rf pulse 
drives them to the 6-state. Consequently, we consider 
two-component bosons trapped in the periodic potential 
formed by interfering laser beams, described by a Bose- 
Hubbard model [lj], 

a={a,b} 

3 \<x,P J 

where c a and S a are the annihilation and creation op- 
erators for states in the internal state a, fi a is the 
chemical potential, Vj j<7 is the external potential with (5, 
the vacuum a-b splitting, absorbed into it, U a p is the 
a state-/3 state on-site interaction strength, and t a is 
the hopping matrix element. The interactions are tun- 



able via Feshbach resonances and spin-dependent lat- 
tices are also available [lOj. For this latter setup, the 
hopping matrix elements may be tuned by the inten- 
sity of the lattices, and introducing small displacements 
of the lattice will reduce the overlap between the Wan- 
nier states of a and b atoms, and therefore may also 
be an efficient way to control the relative size of U aa 
and U a b- The interaction Ubb will be irrelevant: we 
will only consider the case where there is a vanishingly 
small concentration of 6-state particles. In calculating 
the response to RF photons we will take Vj = constant. 
Trap effects will later be included through a local den- 
sity approximation J5|l which is valid for slowly varying 

traps B El EE EE El HE EE EE EE Hp • 

Experimentally the RF spectrum is measured by 
counting the number of atoms transferred from state a to 
b when the system is illuminated by a RF pulse. These 
dynamics are driven by a perturbation 

flrf = ^M/)ri ,,<;,,, • II,-.. (2) 

3 

where j(t) is proportional to the time-dependent ampli- 
tude of the applied RF field multiplied by the dipole ma- 
trix element between states a and b: typically 7 is a 
sinusoidal pulse with frequency uo with a slowly varying 
envelope ensuring a small bandwidth. Due to the small 
wave-number of RF photons, recoil can be neglected. 

For a purely sinusoidal drive, the number of atoms 
transferred per unit time for short times is 

r H = ^E^-^/-^))K/l^fK)l 2 (3) 

where the sum is over the initial states (occupied with 
probability pi = e~@ Ei ) and the final states, all of which 
are eigenstates of H with energies E{ and Ef. We will 
restrict ourselves to T = and the physically relevant 
case where the initial states contain no b- atoms. 



B. Sum Rules 

Taking moments of Eq. [jj 0, EH, the mean ab- 
sorbed photon frequency is 

~JdwY{u) = (4) 

= S-z(t b -t a )f c + (U ab -U aa )g 2 (n). (5) 

We defined 6 to be the vacuum a-b splitting, the local 
phase coherence factor is 

with i and j nearest neighbors, the site filling is n = 
c\c a , and the lattice coordination is z. The zero-distance 
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density-density correlation function is 



92 



C a C a C a c a) 

(n) 2 



(7) 



The second term in Eq. ([5|) may be interpreted as the 
mean shift in the kinetic energy when the spin of an 
atom is flipped. In particular, within a strong-coupling 
mean-field picture (c\ a Cj ja ) = {c\ a )(cj,a) is the conden- 
sate density, which can therefore be measured with this 
technique. The second term in Eq. ([5]) is the shift in the 
interaction energy. 

Our subsequent approximations will satisfy this sum 
rule. This is non-trivial: for example, even in simultane- 
ous limits of = 0, U a b = U aa , and t a — > considered 
in Ref. [6[, their results violate this sum rule by a factor 
of ~ 3. 

Since it plays no role in the remainder of the discussion, 
we will set to zero the vacuum level splitting: S = 0. This 
amounts to working in a "rotating frame" . 



III. RANDOM PHASE APPROXIMATION 



A. General setup and solution 
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FIG. 2: (Color online) Homogeneous system's spectral den- 
sity as a function of cj/U a a and t a /U aa (whiter indicates larger 
spectral density) compared with sum rule prediction (red, sin- 
gle line) . Delta functions are broadened to Lorentzians for vi- 
sualization purposes. (a,b) We take Uba = l.2U aa and tb = t a , 
with (a) ji = 1.98 and (b) ii = 2.02. (c,d) We take parameters 
corresponding to typical 87 Rb experiments: Uba = 1.025U aa 
and h = t a , and take (c) \i = 1.999 and (d) ii = 2.004. In 
both cases, a double peak structure is visible, but the region 
of the phase diagram in which it is important is much smaller 
for 87 Rb parameters than for Fig (a,b)'s parameters. 



To calculate the RF spectrum we employ a time- 
dependent strong-coupling mean-field theory which in- 
cludes k = fluctuations around the static strong- 
coupling Gutzwiller mean field theory 0. This mean 
field theory is exact in the deep Mott limit and in the 
deep superfluid when t a = t&, and it yields fairl> 
rate ground states in the intermediate regime ^ 
[3, 0. Refs. H,[I3| previously used analogous RPA's to 
calculate the Bose-Hubbard model's quasiparticle spec- 
tra and RF spectra with U a b = 0, which reduces to the 
k = single particle spectra. 

We use the homogeneous time-dependent Gutzwiller 
variational ansatz 



W)) = 



(fn(t)\n,0) i +g n (t)\n-l,l) l ) 



(8) 



where \n ai ni ) ) i is the state at site i with n a particles in 
the a state and in the b state. The equation of motion 
for f n (t) and g n (t) are derived by minimizing the action 
S = J dtC, with Lagrangian 



£=(^W-^|ffW-A(#>, 



(9) 



where A is a Lagrange multiplier which enforces conser- 
vation of probablility. At time t = — oo, where j(t) = 0, 
we take g n = 0, and choose f n to minimize (ip\H\ip), 



A/ n = ~t a z (Vnot*fn-i + Vn + laf n+1 ) 
^-n(n- 1) - fin) / n , 



(10) 



where 



a = ^Vnfnfn-i- 



(11) 



Solving the subsequent dynamics to quadratic order in 
7, one finds 



r(f) = N s Jdt'^tMt') X {R) (t-t'), (12) 
where the retarded response function is 

X iR Ht) = -J2^( G n(t)fn-G n (t)f:). (13) 

n 

The Green's functions G n (t) satisfy the equations of mo- 
tion for the g n 's in the absence of an RF field, but in the 
presence of a delta function source, and boundary con- 
dition G n (t) = for t < 0. The relevant equations are 
simplest in Fourier space, where G n (uo) = Jdt e tUJt G n (t) 
obeys 



\fnfn — —LdG n + ^ AnmGn 



(14) 



where A = A + is a Hermitian matrix. The tridagonal 
part A is 



An,n+i = -zt a a^n 

An,n— 1 = 



A r 



zt a a* Vn — 1 

U a 



(15) 
(16) 



- M n-A+^(n-l)(n-2) (17) 
+U ab (n - 1). 



The remaining contribution, 0, is 

@nm = —Ztbfn-lfm-1' 

Specializing to the case where a(t) = ae lUJt , the response 
is given in terms of normalized eigenvectors v m , with 
^ m A nm Vm = CjVn^ • It takes the form of a sum of delta- 
functions, 

7 M = EfevW™^) (19) 

j \ m / 

The / n 's are found at each point in the phase diagram 
by starting with a trial a, solving Eq. (fTUjh then updat- 
ing a via Eq. (jTTJ) and iterating. We find that almost all 
spectral weight typically lies in only one or two peaks. 
Fig. [2] shows sample spectra. The superfluid near the 
Mott state displays a multi-modal spectrum, but in the 
weakly interacting limit only a single peak is seen. An 
avoided crossing is clearly visible in these plots. Fig. [3] 
shows the manifold of spectral peaks in the t a /U aa and 
n/U aa plane, using height to denote frequency and opac- 
ity to denote spectral weight. Taking moments of x R ( UJ ), 
we see that Eq. (|5j) is satisfied. 




FIG. 3: (Color online) Three-dimensional plot of RF spec- 
tral frequencies versus rescaled hopping t a /U aa and rescaled 
chemical potential n/U a a for U a b/U aa = 1.2. Larger opacity 
indicates larger spectral weight. White lines represent con- 
tours of fixed /i, U and oo. The main branch is colored so 
that the progression from green to red to blue corresponds to 
increasing uj. The double peaked spectrum is apparent from 
the "double-valuedness" of the surface. To avoid clutter, nu- 
merical values are omitted from the axes: the Mott plateaus 
occur at frequencies co = 0,0.2U aa and 0AU aa , are each U aa 
wide and the first lobe's critical t is around 0.029£7 a a in 3D. 



B. Limiting Cases 

Although finding the spectrum in Eq. ([T9]) is a trivial 
numerical task, one can gain further insight by consider- 
ing limiting cases. First, when U a b — U aa and t a = tb 
the system possesses an SU(2) symmetry. In this limit 
we find that G n (t) = —iy/nf n Q(t) is constant for t > 0. 
Thus our approximation gives a spectrum which is 
proportional to 5{uo). This result coincides with the ex- 
act behavior of the system: the operator X = . b^ctj is 
a ladder operator, [H : X] = SX, and can only generate 
excitations with energy S (set equal to zero in our calcula- 
tion) . The fact that our approximations correctly capture 
this behavior is nontrivial: in a field theoretic language 
one would say that our equation of motion approach in- 
cludes the vertex corrections necessary for satisfying the 
relevant "Ward identities" [II [II, [H . 

The current 87 Rb experiments are slightly perturbed 
from this limit, with r] = (U a b — U aa )/U aa ~ —0.025 
and tb = t a . We find that the J-function is shifted by 
a frequency proportional to ry, but that the total spec- 
tral weight remains concentrated on that one frequency: 
the sum of the spectral weights at all other frequencies 
scale as rj 2 . Consequently it is an excellent approxima- 
tion to treat the spectrum as a delta-function, and our 
RPA calculation reduces to the results in [5[ . We empha- 
size however that other atoms, such as Cesium, can be 
in a regime where r] is large. 

We gain further insight by considering the superfluid 
near the Mott phase with t a /U a <C 1. Here one can trun- 
cate the basis to two states with total particle number n 
and n + 1 on each site. Then the / n 's and G n 's can be 
found analytically: one only needs to solve 2x2 linear 
algebra problems. In the tb = 0, U a b = U aa limit, this 
is similar to Ref. @'s approach, but includes the hop- 
ping self consistently, allowing us to satisfy the sum rule 
Eq. ([5]). This truncation is exact in the small t a limit, 
and yields 

X (jR) M = A + 5(uj-uj + )+A_8(uj-uj-) (20) 

with 

u± = ^±^A 2 + (^) 2 (21) 

where 

Cl = (Uab-Uaa)(n-l)-\-Zt a fn + l(n+l) 

e 2 = (U ah -U aa )n + z[t a (n + l)-t h )fl 

A = -Vn(n+l)t a zf n f n+1 . (22) 

if n > 1 and 

ei zt a fi 

£2 = z(t a - t b )fo 

A = 0. (23) 

if n = (here, only the 62 peak has non-zero spectral 
weight). We omit the cumbersome analytic expressions 
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for the spectral weights A±. The spectrum consist of 
two peaks - hybridized versions of the excitations car- 
icatured in Fig. [TJ One can identify e\ and 62 as the 
energies of those caricature processes, recognizing that 
the hybridization term, A, grows with t a . The avoided 
crossing between these modes is evident in Fig. El 



C. Inhomogeneous spectrum 

We model the trapped spectrum through a local den- 
sity approximation. We assume that a given point in the 
trap has the properties of a homogeneous gas with chem- 
ical potential fi(r) = [Iq — V(r). In Fig. [H we show the 
density profile and the spectrum corresponding to each 
point in space. Also shown is the trap averaged spec- 
trum. The bimodality of the homogeneous spectrum is 
quite effectively washed out by the inhomogeneous broad- 
ening of the trap. On the other hand, if one images the 
atoms flipped into the b state as in Ref. [3], there is a 
clear qualitative signature of the bimodality. If one ex- 
cites the system with an RF pulse whose frequency lies 
between the resonant frequencies of two Mott plateaus, 
one will excite two "shells" of atoms. These shells should 
be clearly visible, even in column integrated data. 



IV. CONCLUSIONS AND DISCUSSION 

In this paper we have shown that the RF spectra of 
a homogeneous Bose gas in an optical lattice will have 
two (or more) peaks in the superfluid state when the 
parameters are tuned close to the superfluid-Mott insu- 
lator phase transition. Physically, this bimodality is a 
result of the strong correlations in the system. These 
correlations result in two distinct forms of excitations 
(which are strongly hybridized): those involving "core" 
atoms, and those involving delocalized atoms. When 
77 = (U a b — U aa )/U aa is small, such as in the experiments 
on 8T Rb, this bimodality is absent. 

Our approach, based upon applying linear response to 
a time dependent Gutzwiller mean field theory, is both 
simple and quite general. It allows arbitrary interactions 
between both spin states, and it allows arbitrary spin- 
dependent hopping rates. The major weakness of the 
theory is that it fails to fully account for short range- 
correlations: the atoms are in a quantum superposition 
of being completely delocalized, and being confined to 
a single site. The physical significance of this approxi- 
mation is most clearly seen when one considers the case 
where the final- state atoms have no interactions, U a b — 0, 
and see no trap or lattice. Imaging the 6-atoms after a 
time-of-flight is analogous to momentum resolved pho- 
toemission [25| , and would reveal the dispersion relation- 
ship of the single-particle excitations. The fact that the 
spectrum consists of two sharp peaks means that all of 
the non-condensed atoms are approximated to have the 
same energy. One will also see that their momentum is 
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FIG. 4: (Color online) (a) Density n as a function of dis- 
tance to trap center rescaled by the lattice spacing, r/d, in 
a local density approximation. For all subfigures, we take 
ta/Uaa = 0.004, which is moderately smaller than the tip of 
the first Mott lobe, (b-e) Left: spectrum of a homogeneous 
gas with density n(r), representing the spatially resolved spec- 
trum observed in an experiment on a trapped gas. Horizontal 
axis is position, vertical is frequency, color from dark to light 
represents increasing spectral density. Continuous (red) curve 
denotes sum rule result for (u). We round the ^-functions to 
Lorentzians for visualization. Right: trap- averaged spectrum 
for a 3D trap within our RPA (black, solid line) compared 
with sum rule (red, dashed line), (b) U a b — 1.2U aa ,tb = t a (c) 
Uab = U aa ,t b = t a +0.1U aa (d) U ab = 1.2U aa ,t b = t a + 0.1U aa 
(e) 87 Rb parameters: U a b = l.Q25U aa ,tb = ta- 
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uniformly distributed throughout the first Brillioun zone. 
In the strong lattice limit, where the bandwidth is small, 
this approximation is not severe. 
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